function Mesh = process_mesh_2d(PDE, h, basis_type)
Mesh.basis_type = basis_type;
Mesh.left = PDE.domain(1);
Mesh.right = PDE.domain(2);
Mesh.bottom = PDE.domain(3);
Mesh.top = PDE.domain(4);
Mesh.Nx = (Mesh.right-Mesh.left)/h;
Mesh.Ny = (Mesh.top-Mesh.bottom)/h;
[node, elem, Mesh.bde] = process_square_2d(Mesh.left, Mesh.right, Mesh.bottom, Mesh.top, Mesh.Nx, Mesh.Ny);
Mesh.Np = (Mesh.Nx+1)*(Mesh.Ny+1);
[Mesh.node, Mesh.elem, Mesh.bdn] = process_grid_interpolation(node, elem, Mesh.bde, basis_type);
Mesh.Ne = size(Mesh.elem,1);
Mesh.Nb = size(Mesh.node,1);
Mesh.Nlb = size(Mesh.elem,2);
end

%% process_square_2d
function [node, elem, bde] = process_square_2d(left, right, bottom, top, Nx, Ny)
[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
Ne = size(elem,1);
bd_bottom = [1:2*Ny:Ne;
    1:Ny+1:Np-Ny-1;
    Ny+2:Ny+1:Np]';
bd_right = [(Nx-1)*2*Ny+2:2:Ne;
    Np-Ny:Np-1;
    Np-Ny+1:Np]';
bd_top = [2*Ny:2*Ny:Ne;
    2*Ny+2:Ny+1:Np;
    Ny+1:Ny+1:Np-1]';
bd_top = sortrows(bd_top,1,"descend");
bd_left = [1:2:2*Ny-1;
    2:Ny+1;
    1:Ny]';
bd_left = sortrows(bd_left,1,"descend");
bde = [bd_bottom; bd_right; bd_top; bd_left];
end

%% boundaryedge
% function bde = boundaryedges(elem, bdary)
% E_circle = [elem, elem(:,1)];
% edges = [E_circle(:,1:2); E_circle(:,2:3); E_circle(:,3:4)];
% [~, loc] = ismember(edges, bdary, 'rows');
% E_bd = reshape(loc, size(elem,1), 3);
% [rows, cols] = find(E_bd);
% index0 = zeros(size(rows,1), 1);
% for i = 1:size(rows,1)
%     index0(E_bd(rows(i), cols(i))) = rows(i);
% end
% bde = [index0, bdary];
% end

%% process_grid_interpolation
function [new_node, new_elem, new_bdn] = process_grid_interpolation(node, elem, bde, basis_type)
bdn = unique(bde(:,2:3), "stable");
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
        new_bdn = bdn;
    case "P1b"
        P = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            P(k,:) = polygon_centroid(V);
        end
        new_col = (size(node,1)+1:size(node,1)+size(elem,1))';
        new_node = [node; P];
        new_elem = [elem, new_col];
        new_bdn = bdn;
    case "P2"
        v0 = reshape(elem', 3*size(elem,1), 1);
        v1 = reshape(([elem(:,2:3), elem(:,1)])', 3*size(elem,1), 1);
        totalEdge = [v0,v1];
        index1 = zeros(size(totalEdge,1),2);
        k = 1;
        for i = 1:size(totalEdge,1)
            if ~all(totalEdge(i,:) == [0, 0])
                for j = i+1:size(totalEdge,1)
                    if all(totalEdge(i,:) == flip(totalEdge(j,:)))
                        totalEdge(j,:) = [0, 0];
                        index1(k,:) = [j, i];
                        k = k+1;
                    end
                end
            end
        end
        uniqueEdge = totalEdge(~(totalEdge(:,1) == 0),:);
        P = (node(uniqueEdge(:,1),:) + node(uniqueEdge(:,2),:))/2;
        index0 = setdiff((1:size(totalEdge,1))', sort(index1(:)));
        index1 = sortrows(index1(~(index1(:,1) == 0),:));
        index2 = [find(totalEdge(:,1) ~= 0), (1:size(P,1))'];
        new_col = zeros(size(elem,1), 3);
        k = 1;
        for i = 1:3:size(totalEdge,1)-2
            new_col(k,:) = [i, i+1, i+2];
            k = k+1;
        end
        for k = 1:size(index1,1)
            idx = (new_col == index1(k,1));
            new_col(idx) = index1(k,2);
        end
        for k = 1:size(index2,1)
            idx = (new_col == index2(k,1));
            new_col(idx) = index2(k,2);
            idx = (index0 == index2(k,1));
            index0(idx) = index2(k,2);
        end
        [~,index3] = sortrows(bde(:,1));
        [~,index3] = sortrows(index3);
        index0 = index0(index3);
        index0 = index0 + size(node,1);
        new_node = [node; P];
        new_elem = [elem, new_col+size(node,1)];
        new_bdn = reshape([bdn'; index0'],[],1);
    otherwise
        error("Invalid basis type.");
end
end

%% polygon_centroid
function center = polygon_centroid(V)
n = size(V,1);
A = 0;
center_x = 0;
center_y = 0;
for i = 1:n
    j = mod(i,n) + 1;
    xi = V(i,1);
    yi = V(i,2);
    xj = V(j,1);
    yj = V(j,2);
    factor = xi * yj - xj * yi;
    A = A + factor;
    center_x = center_x + (xi + xj) * factor;
    center_y = center_y + (yi + yj) * factor;
end
center_x = center_x / (3 * A);
center_y = center_y / (3 * A);
center = [center_x, center_y];
end